Liquid Conservation and Non-local Interface Dynamics in Imbibition 
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The propagation and roughening of a liquid-gas interface moving through a disordered medium 
under the influence of capillary forces is considered. The system is described by a phase-field model 
with conserved dynamics and spatial disorder is introduced through a quenched random field. Liquid 
conservation leads to slowing down of the average interface position H and imposes an intrinsic 
correlation length £ x ~ H 1 ^ 2 on the spatial fluctuations of the interface. The interface is statistically 
self affine in space, with global roughness exponent \— 1-25 and exhibits anomalous scaling. 



;.35.Ct, 



.35.Fx, 47.55.Mh, 05.70.Ln 



The dynamics of an invading liquid front in a disor- 
dered medium involves many problems of interest in sta- 
tistical physics such as transport and diffusion in random 
media [EU3], roughening of a moving interface j^], and 
pinning [0] if the interface stops. Many experiments have 
been conducted on such phenomena , examples include 
Hele-Shaw cells B and the spontaneous imbibition of wa- 
ter in paper |7| pd|. In the latter case, capillary forces 
drive the liquid until balanced by loss of water by evapo- 
ration or hydrostatic pressure. Often local theories such 
as directed percolation depinning have been success- 
ful in describing the properties of such pinned fronts . 

Local interface theories cannot however provide a com- 
plete description of imbibition since the transport of liq- 
uid to the front from the reservoir must be taken into 
account. In particular, this leads to a slowing down of 
the invading front even without evaporation or gravity 
|^2| . In this letter a model of imbibition that explic- 
itly addresses this issue through the inclusion of liquid 
conservation, is introduced. This is achieved by a gen- 
eralised Cahn-Hilliard equation with a non-equilibrium 
boundary condition that couples the system to a liquid 
reservoir. Two different cases are considered: a liquid 
front invading a disordered medium without gravity or 
evaporation, and a steady state front maintained at a 
constant distance from the liquid reservoir as done in the 
experiments of Horvath and Stanley M] . 

An obvious consequence of liquid conservation is that 
the interface dynamics is nonlocal. One manifestation 
of this non-locality is that the average interface height 
H grows as i 1 / 2 as might be expected by Washburn's 
equation Jl2] (the velocity, H ~ 1/H). The interface is 
found to be superrough, with global roughness exponent 
X — 1.25, and exhibits anomalous scaling |} 13|. Interface 
fluctuations saturate at a correlation length £ x , which 
is controlled by the average interface height such that 
£ x ^H 1 / 2 . This implies that the average interface width 
Hi- The relationship between £ x and H is a di- 



rect consequence of the conservation law and arises from 
the interplay between driving force and curvature. No 
such relation occurs in local models S. The temporal 
correlations of the interface qualitatively agree with the 
experiments of Ref . pj . It is further found that the front 
moves by avalanches, as reflected in the behaviour of the 
higher moments of the temporal correlations (r| . 

Phase field model of imbibition: To describe the liq- 
uid/gas system a locally conserved field <f>(x,t) is intro- 
duced such that <f> = +1 and —1 are the liquid and gas 
phases respectively. The free energy functional of the sys- 
tem can then be written, F{<j>} = / dx[(V(j)) 2 /2 + V(0)], 
where V(x, <j>) = -0(x, t) 2 /2 + 0(x, f) 4 /4 - a(x)0(x, t). 
The liquid/gas interface is thus sustained in the stan- 
dard manner by the gradient energy term and double 
well structure of V. The variable a(x) accounts for 
the random nature of the medium and has correlations; 
(a(x))=a, and (a(x)a(x')) -a 2 = (Aa) 2 <5(x-x'). The 
propensity of the medium to absorb liquid can be con- 
trolled by a since it determines the local chemical po- 
tential (/i = —5F/6(f>). The dynamical evolution of <f> is 
determined by a continuity equation 9t<^+V-j = 0, where 
j(x, t) = — V/x(x, t). The resulting equation of motion, 
i.e., 
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(x)] 
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is a slightly modified Cahn-Hilliard equation |16fl . The 
absorption of liquid by a dry random medium can be 
modeled with the appropriate boundary conditions. To 
reflect the presence of an infinite resevior at y < the 
chemical potential at y = is fixed to be zero. To induce 
the propagation of the front, the top end of the system is 
kept "dry" (i.e., <fi(y — > 00) = — 1) and a > 0. Equation 
(|l|) can be augmented to include evaporation by adding 
a term proportional to (0+1). 

This model was specifically constructed to incorporate 
several essential physical features of imbibition. More 
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precisely the model incorporates the conserved nature of 
the fluid which leads directly to a Gibbs-Thomson bound- 
ary condition on the liquid front and an average front 
velocity that decreases in time [|l7| . To see these proper- 
ties and gain more insight into the model it is useful to 
project out the equation of motion of the interface using 
standard projection operator methods ]l8| , ^9| . Expand- 
ing to lowest order in the interfacial curvature (/C) and a 



gives; 



Acj) 



dx'Q(x, h\x' , h')h' = 7](x, h) + oK 



(2) 



where Q is the half-plane Green's function given by 



G(x,y\x',y') 



(x-x') 2 + (y-y') 2 



which reflects a 



(x-x') 2 + (y+y') 

broken translational symmetry along the y-axis due to 
the reservoir. The interface is defined by the single- 
valued function h{x,t), a = J dy(<pQ) 2 is the surface ten- 
sion, 0o is the one-dimensional kink solution of Eq. (|l|) 
and A(j) = 2 is the miscibility gap. The quenched noise 
can be written rj{x,h) = J dy<p' (y — h(x, t)) a(x, y) ~ 
2a{x, h) in the sharp interface limit. Similar non-local 
equations also arise in the context of Laplacian fluid flow 
|plf and step growth 0). Previous theoretical work on 
equations for interfaces at the depinning transition 
and nonlinear equations with long range kernels 23 does 
not apply to the situation encountered here. 

The Gibbs-Thomson condition, n\i n t ~ IC can be im- 
mediately obtained from Eq. (^) in the limit h = since 
r\ is the chemical potential at the interface. In addition 
since is a conserved field it is easy to show that the 
normal interface velocity is proportional to the gradient 
of /i at the interface, i.e., Atfiv n — —d n ^i\i n t- Further in- 
formation can be obtained by linearising in h to obtain; 



Lateral correlation length: Another useful result com- 
ing from Eq. (^) is the existence of a lateral length scale 
separating two different modes of damping of the in- 
terface fluctuations. Long wavelength fluctuations are 
damped by the advancing motion of the front (i.e., 
\k\Hhk), while short wavelength fluctuations are damped 
by surface tension (i.e., a\k\ 3 hk)- The length scale sep- 
arating these two modes £ x ~ (er/ij) 1 / 2 = (aH/a) 1 / 2 is 
closely related to the Mullins-Sekerka instability of driven 
Laplacian fronts p5[ , although the situation is reversed 
here. Because fluid is transported towards the front from 
behind, advanced (retarded) parts receive less (more) 
mass than the average and the front is stabilised at long 
length scales. 

This result can be understood in more physical terms 
as follows. Due to the Gibbs-Thomson effect a local 
"bulge" of vertical extent W and lateral size £ alters the 
chemical potential by A/i ~ aW/^ 2 . On the other hand, 
the average gradient in /i in the bulk liquid induces a dif- 
ference A/i ~ aW/H across a vertical distance W. These 
two differences balance each other at a length given by 



h k (l 



-2\k\H 



where h k are the Fourier components of h and H = Iiq is 
the average interface position J|4| . Equation (^) can now 
be used to obtain an expression for the average interface 
position (Washburn's equation) and lateral correlation 
length as follows. 

Washburn's equation: In the limit k — > 0, Eq. (||) re- 
duces to H — a/(2H). Thus a planar interface, without 
evaporation is never pinned and moves as H (t) — (at) 1 / 2 . 
This is a standard result known as Washburn's equation 
fl2f and can be understood as follows. In the absence of 
a liquid reservoir \x = —a throughout the system. When a 
reservoir is included at the lower boundary fj, is set to zero 
at y = 0. For a slowly moving front /i(y) satisfies a dif- 
fusion equation in the bulk phases with quasi-stationary 
solution fi w —ay/H for y < H, and /j,(y) = —a for 
y > H. This result is equivalent to Washburn's equation 
since H = v n = — (9 n /i|i„t)/A0. In essence the gradi- 
ent in /i creates a current from the reservoir towards the 
interface, causing it to advance. 



(4) 



By this mechanism, one expects the front to be smoother 
on length scales larger than £ x as compared to smaller- 
scales. While these interface results provide useful insight 
into the front propagation the behavior was examined in 
more detail by direct numerical simulation. 



) + \k\H h k (l + e - 2 l fc l«) = \k\ ( Vk - ak 2 h k ) 

(3) 
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FIG. 1. Configurations of the rising interface at time inter- 
vals At = 10 3 . For this system, the average height H ~ t 0A9 , 
and the width W(t) ~ t a32 , according to least-square fits. 
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Numerical integrations: Equation (|l]) was integrated 
on a square grid with spacing Ax = 1 using various lat- 
eral system sizes L x . The disorder a(x) was a random 
variable on each grid point. Different distributions (ex- 
ponential, Gaussian, etc.) gave the same results. The 
interface position h(x, t) was determined by a linear in- 
terpolation of the zero of the phase field, 4>(x, h(x, t)) = 0, 
with overhangs (hardly present here) cut off. 

Fig. |l| shows the rise of an initially flat front. Fronts 
are shown at equal time intervals, displaying Washburn 
behavior H 2 (t) = at, also shown in the Figure. The 
vertical and lateral extent of interface fluctuations in- 
crease with height. The early time behavior of the width 
W(t) = {(h(x,t) — h) 2 ) 1 / 2 , where overbar denotes a sys- 
tem average and the brackets average over the disorder, 
is also shown in Fig. [|. It is found that W(t) ~ t® with 
(3 = 0.32 ±0.02. 

In a recent experiment of paper wetting JhJ , a steady 
state situation was reached by pulling the paper sheet at 
a constant velocity v = — vy towards the reservoir. To 
model this case, Eq. (w) can be written as 



d t 



-V 2 [V 2 



a(x — vi) 



(5) 



keeping the original boundary conditions. In this case, 
the interface remains at a fixed average height H — a/2v, 
at which a freely rising interface would have velocity v. 



(a) 




(b) 




FIG. 2. (a) Structure factors for systems of height 
# = 50, 100, 150, 200 and L x — 2H. fb) Correlation functions 
G2(r,H), scaled according to Eq. M) for various v ~ H" 1 ^ 2 
and Aa. (c) Comparison of the correlation functions in the 
steady state, G2(r,H) (solid lines) for i/ = 25,50, 100 and of 
the rising case G2(r,tii) (dashed lines), at times tn — H 2 /a. 

Equation (^|) was also integrated numerically. The 
role of £ x as a maximal range of correlated roughness 
becomes apparent in the steady-state structure factor 
S(k,H) = {\h(k,t)\ 2 ) and spatial height difference cor- 
relation function G 2 (r,H) = {(h(x + r, t) - h(x, t)) 2 ) 1/2 
p6[ , both shown on Fig. |2| The lateral fluctuations of the 
interface are correlated only up to £ x < L. The global 



roughness exponent x is taken from the power law de- 
cay of S(k > £x\iT) ~ fc-( 2x+1 ) = fc~ 3 - 5±0 - 2 , yielding 
X ~ 1.25. The relation £x ~ u~ 1/2 ~ (H/a) 1 / 2 is con- 
firmed in Fig. ^, were the spatial correlation function is 
rescaled as 



G 2 (r, H) = Aa «-* /2 g (r v 1/2 ) 



(6) 



for different v (i.e., different H) and Aa. The scaling 
function g(u) is constant for b» 1 and g{u) ~ u Xloc if 
m«1, with xioc — 1- The interface is thus superrough, 
and G2 (r, H) shows anomalous scaling in the sense that 
Xioc<X- G2{r,H) increases with H for all r p|Jl3[]. 

Quasi- stationary rise: The numerical data indicate 
that spatial correlations in the interface position are the 
same for both sets of simulations. This can be seen in 
Fig. H where steady state correlation functions at fixed 
heights H are compared with the correlation functions 
of a freely rising front, obtained from Eq. (|l|), at times 
tn = H 2 1 'a. This implies that the interface is always in a 
quasi-stationary state in the freely rising case. The scal- 
ing form, Eq. (||), is thus also valid for the freely rising 
front provided that the time dependence H(t) = (at) 1 / 2 is 
used, such that £ x ~ (t/a) 1 ^ |19|. This indicates that the 
dynamical correlation length ^f~t 1/,z increases at a rate 
faster or equal to t 1 / 4 , so that the interface fluctuations 
can always catch up with the maximal correlation range. 
This also implies W~£* -i x/4 = f 3 , with ~ 0.31, in 
agreement with Wit) shown in Fig. |l|. 
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FIG. 3. Steady state temporal correlation functions for 
systems of size L x = 400, at heights H = 50, 100, 200. In 
the inset, the moments q — 2, 4, 6 (from bottom to top) for 
H = 100 are compared. C'2 has exponent (3% = 0.85, while for 
the higher moments, the effective exponents /?4 ~ 0.76 and 
fa » 0.69. 

Temporal correlations: To compare with the ex- 
perimental work of Ref. the temporal cor- 
relations in the steady state, obtained from Eq. 
(||), were analysed through the function C q (t,H) = 
(\h(x,t + s) — h(x,s)\ q ) 1 / q also taking higher moments 



3 



q > 2 into account. The early time logarithmic slopes 
of the higher moments (/3 q ) decrease with q, as shown in 
the inset of Fig. |[ in contrast to the spatial G q for which 
X q is independent of q. This lack of scaling reflects the 
propagation of the front by avalanches (easily visible in 
Fig. [I]), also found for local interfaces at the depinning 
transition p5| , pj| . 

Nevertheless, our model agrees with the experiments 
in three ways, visible in the main part of Fig. |^: (i) 
the P q are independent of H; (ii) the late time satura- 
tion level of C q increases with H indicating larger overall 
roughness, and (in) C q (t,H) at fixed small t decreases 
with H , indicating faster intrinsic avalanche velocity for 
smaller H . Quantitatively our result /?2 = 0.85 ± 0.04 
differs from the experimental value 0.56, but so does the 
average interface velocity v ~ H~ l e in Ref. [fl0[ com- 
pared to Washburn behavior in our model [B7L Since 
water flow dynamics in real paper can be complicated 
by various phenomena such as e.g. fiber swelling the 
disagreement is not too surprising and should be clarified 
in future investigations. 

Conclusion: We have constructed a simple but flexi- 
ble model for spontaneous imbibition of a liquid from a 
reservoir into a disordered medium. Conservation of liq- 
uid leads to Washburn's relation for the average height 
H ^i 1 / 2 [jl2), a behavior which does not occur naturally 
in local interface models. The interface dynamics be- 
comes nonlocal, and an intrinsic length scale £ x emerges 
as given by Eq. ([!]). This length scale is a distinct predic- 
tion of the model, and intimately relates the fluctuations 
of the interface with its average progression. The ex- 
istence and dynamical increase of this region should be 
experimentally easier to observe than precise values of 
scaling exponents on scales r < £ x . From the model, 
spatial scaling is found to be anomalous with a global 
roughness exponent \ — 1-25, similar to what has been 
seen in other models of driven interfaces in disordered 
media [^8| . The growth exponent (3 is explained in terms 
of £ x and \- Acknowledgements: This work has in part 
been supported by the Academy of Finland and a Re- 
search Corporation Grant, No. CC4787 (KRE). 
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